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SUMMARY 

This report describes the methods used to compute the mass, struc- 
tural stiffness and aero<fynamic forces in the form of influence coefficient 
matrices as applied to a flutter analysis of the DAST ARW-1 wing. The 
DAST wing was chosen since wind tunnel flutter test data and zero speed 
vibration data of the modes and frequencies exist and are available for 
comparison. 

The report also contains a derivation of the equations of motion that 
can be used to apply the modal method for flutter suppression. A compari- 
son of the open loop flutter prediction with both wind tunnel data and other 
analytical methods is presented. 
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INTRODUCTION 


Real time, feedback control for flutter suppression is under serious 
study and consideration for aircraft (Refs. 1,2, 3, 4). The modal method 
(Refs. 5,6) is well suited for application in flutter suppression since the 
onset of flutter may be adequately described as a linear system instability 
(Ref. 7). Previous analyses of the problem have employed generalized 
coordinates, based on zero airspeed vibration modes or other fixed wing 
deformation shapes, from which generalized aerodynamic forces have been 
computed (Refs. 1,2,4). The contribution of this report is that physical 
coordinates of bending and torsion of the wing structure are directly employed, 
and that constant influence coefficient matrices are used to describe the struc- 
tural, inertial and aerodynamic forces over a wide rar^e of Mach numbers and 
airspeed. The aerodynamic influence coefficients were obtained through a 
modification of the SOUSSA digital program (Refs. 8,9) generated with the 
assistance of Prof. L. Morino of Boston University, Dr. E. C. Yates and 
H. Cunningham of NASA/LaRC. In contrast, the aerodjmamic coefficients used 
in Ref. 1,4 were obtained using a doublet lattice method (Ref. 10). The method 
of Pade^ approximants (Refs. 7, 11, 12, 13) was applied to derive the aerodynamic 
influence coefficients in the real time domain. The structural influence coef- 
ficients were obtained through the use of the SPAR computer program at 
NASA/LaRC. Finally, the structural and geometric data of the DAST ARW-1 
Wir^, at 111 grid points, was supplied by Mr. R. Doggett of NASA/LaRC. 


This study was funded by NASA Langley Research Center under Contract 
NASI- 15593. 
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I. Description of the Wing Flutter Model 

The DAST ARW-1 wing was designed at the Langley Research Center as 
a swept back, cantilevered wind tunnel flutter model of a protot3rpe» remotely- 
piloted, drone aircraft used to study active control concepts including flutter 
suppression. The data used in this report for determining the inertial and 
structural characteristics of the wing, as well as the results of vibration and 
wind tunnel flutter tests, was furnished by Mr. R. Doggett of NASA/LaRC. 

The geometric planform and dimensions of the wing are shown in Fig. 1. 
The leading edge has a sweep back angle of 44. 32°. The wing has a taper ratio 
of . 392 and an aspect ratio of 6. 4. The airfoil is a NACA 65A10 section. The 
main structural beam is a single tapered aluminum bar construction with a 
cruciform cross section (see Fig. 2). The dimensions of the spar cross section 

at various locations along the length are shown in Table 1. The measured stiff- 
ness distribution is shown in Fig. 3 in terms of the bending and torsional 
stiffness, El and GJ curves. 

The wing is divided into eight pod sections by means of seven ribs oriented 
in the stream direction (see Fig. 4). Each section contains concentrated masses 
rigidly connected to the main beam to provide realistic mass offsets with re- 
spect to the local elastic axis. Each section is covered with balsa inserts and 
the aerodynamic shape is maintained by a precured fiberglass cover. 

A control surface is provided along the trailing e(fee, equipped with an 
electro-hydraulic servo-actuator. The surface hinge line is located at 80% of 
the local streamwise chord. The reaction torques of the actuator are con- 
strained by a link to the main structural beam in the control surface pod section. 

A Cartesian coordinate system is used in the analysis. The origin of the 
system is at the intersection of the wing root chord and the wing leading edge. 
The x-axis is positive forward in the streamwise direction. The z-axis is 
positive down, and the y-axis forms a right hand system. (See Fig. 4). The 
structural axis of the beam defines the y* coordinate, rotated at an angle of 
40.7° with respect to the y axis. The origin of the x' , y' , z’ system is at 
(-.372364, 0., 0.). 
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Thus 



where 6 = 40.7°. 


The dynamical coordinates are located alor^ the y* axis (see Fig. 5). 
These consist of seven vertical deflections, h (y’p (i = 1,7), seven rotations 
in the air stream direction, a (y'.) (i = 1,7), and a single rotation about the 
control surface hinge line, 6 . Thus the dynamical coordinates form a 
15 X 1 vector given by 


r 



C w } = 



( 2 ) 




»7 

6 .y 


The equations of motion will be generated in terms of the forces and moments 
affecting these 15 degrees-of-freedom. The center line, or root section, is 
constrained to maintain zero deflection in the vertical direction and in stream- 
wise rotation. Thus h (0) = a(0) = 0. 
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II. structural Influence Coefficients 


The structural influence coefficients were computed utilizing the SPAR 
computer program at the LRC computer facility. The SPAR program re- 
quires that the wing be decomposed into a series of grid points. At each 
point, six degrees-of-freedom are permitted. There are three translations 
along, and three rotations about each of the three axes. A finite element 
method is employed to compute the linear relationships between the defor- 
mations of the grid points with respect to one another and the resulting forces 
and moments resisting these deformations. For the DAST wing, a grid 
decomposition of 111 was used, resulting in 666 degrees-of-freedom. (See 
Figure 6). Of these, points 103 through 109 correspond to the concentrated 
masses which undergo rigid motions without relative structural deformations. 
Point 111 is the control surface linkage constraint for rotation about the hinge line. 

To produce the structural influence coefficients for the 15 degrees-of- 
freedom defined by the (fynamical coordinate vector, w , we constrained 7 
points in 2 degrees-of-freedom (vertical deflection and rotation in the flight 
direction) and the 15th degree-of-freedom to be a pure rotation about the 
hinge line. All other degrees-of-freedom In the original SPAR deformation 
are left unconstrained except for the root section grid points (97, 100 and 102) 
which are constrained in three directions. (See Figure 7). The SPAR program 
solves a static equilibrium problem for which one coordinate of w is set equal 
to unity, and the other twenty-three are constrained to be zero. The forces 
and moments at the 24 locations are computed by solving a set of 24 equations 
of equilibrium. Thus, we have 




M. 
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PC 1c *k 
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where the -tth element of w , w , = 1 and all other elements, including 
h(0) and a(0), are set equal to zero. The forces, F . , and the moments, 

U j 'v 

M , are the reactions at the root section required to hold the root section 

U j V 

undeformed. 

The elements of the -tth column of the influence coefficient matrix (K) 
are given by 



The deflections at the remaining 642 grid points are left unconstrained. 
The matrix of influence coefficients is given in Table 2. The units are in 
Newtons per meter of deflection. Newtons per radian. Newton-meters per 
meter of deflection and Newton- meters per radian arranged as follows : 


I 

Newtons/meter (7x7) | Newtons/rad (7x8) 


I 

Newton-meter/meter j 
(8x7) I 


Newton-meter /rad 
(8x8) , 


(4) 
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ni. Mass Data 


The SPAR program carries out a vibration analysis of the DAST wing 
by solving an eigenvalue problem utilizing the input mass and stiffness data at 
the 111 grid points. Of the 666 degrees-of-freedom, 18 at the root section 
(grid points 97, 100 and 102) are constrained to be fixed, and 42 correspond 
to the 7 rigidly attached concentrated masses (103 through 109). For the 
purposes of this study, the number of degrees-of-freedom has been reduced 
to 15 plus 9 fixed, root degrees-of-freedom. In order to produce a simulation 
in which the significant vibration modes are well represented, it is necessary 
to compute a set of lumped masses at the c.g. *s of the seven sections, which, 
together with the rigid concentrated masses and the 15x15 influence coefficient 
matrix, will reproduce the significant low-order vibration modes. 

To accomplish this, the distributed beam and plate masses have been 
summed in 7 sections and are listed in Table 3. The concentrated masses are 
listed in Table 4. 

The combined masses acting at each of the sections’ c.g. 's form the 
di£^onal elements of the mass matrix (in units of kilograms in the mks 
system). The off-diagonal elements are obtained from the static unbalance 
due to the offset of the concentrated masses from the coordinate c.g. ’s (see 
Table 4). M , the desired mass matrix used in the equations of motion, is 
listed in Table 5. 
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IV. Frequencies and Modes 


To ensure that the SPAR representation of the DAST cantilevered 
wing is a valid simulation of the wind tunnel model, a comparison of the first 
four computed and measured vibration modes was carried out using the full 
111 grid-point model. Figure 8 contains the results of the vibration test fur- 
nished by R. Doggett of NASA/LaRC. Figures (9a), (9b), (9c), and (9d) are 
plots of the same modes with the SPAR program using all of the 111 grid 
points. The comparison is seen to be good. Figures 10 and 11 contain the 
next two highest modes obtained by the SPAR program for which no vibration 
data is available. 

Finally, in order to test the validity of the reduced 15 degrees-of- 
freedom model, an eigenvalue analysis using the mass, M , and influence 
coefficient, K , matrices was carried out. The results of the analysis are 
shown in Figures (12a), (12b), (12c), and (12d) for the first four modes. The 
agreement is seen to be good. The 15 degrees-of-freedom eigenvalues are 
shown in Table 6 together with the vibration test frequencies and those fre- 
quencies obtained with the full 666 degrees-of-freedom. 
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V. Equations of Motion 

The forces acting on the wing in the air speed region containing the 
flutter speed are assumed to consist of the followii^: 

(a) inertia 

(b) unsteady aerod 3 mamic forces 

(c) structural restraint to wing deformation 

(d) random aerodynamic forces due to wind gusts 

(e) a stabilizing feedback torque acting on the control response 

(f) aerodynamic forces due to the control surface deflection 


To simulate these forces, assuming small deflections, we require the dynamic 
coordinate vector, w , its first and second time derivatives, w and w , the 


unsteady lift and moment vector, x_ , its time derivative vector, x , 

^15x1 P 

three scalar gust variables and w^), and the scalar control torque. 


The equations of motion are given by 
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where 


I 

D 

w 

M 

K 


H 


1 


f 




B 


15x15 identity matrix 

15x15 structural dampii^ matrix (see Table 8) 

15x15 matrix of masses 

15x15 structural influence coefficient matrix from the 
SPAR program 

H are constant 15x15 aerodynamic influence coefficient 

O 

matrices obtained by Fade Approximates from the SOUSSA 
output. (See Table 7b, c,d) 

15x1 vector of O's except at the points of application of the 
feedback torque, u^ , where aJid b^^^ = 1. 


ro 

0 

0 

0 

0 

0 

0 



15x1 


(6b) 


H 
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c 

q 

V 

p 

b 

w 

g 


15x15 steady state lift and moment distribution matrix 

obtained from SOUSSA (See Table 7a) 

half the mean aerodynamic chord (. 2524379 m) 

, 2 

2 p V (dynamic pressure) 

air speed 

density 

reference length used in the SOUSSA program (b = 1 inch) 
scalar random wind gust 
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To facilitate the transformation to modal coordinates, the matrix con- 


taining M and H may be inverted. The, Equation 6 becomes 

O 


r ^ 

w 


-t < 


w \ = 


X 



-M~^D 


^ -1 
1. M D 
3 w 




(6c) 


r 


r ^ 

0 



-1 

>“a M 

, -1 

/MB 
1 -1 1 

M y 
-1 

Hg M b’ 


Hg M y 





(.13-^qXjjj + .565-2-x2p) 
c 


To obtain a good model for the structural damping matrix, , we 

make use of the approximation that the structure provides approximately . 5% 
of critical damping. 


Thus, let 


-1 2 -1 
M K = U ( 03. ) U 


(6d) 


where the matrix U is the matrix of eigenvectors of M K and 03. is 
the diagonal matrix of the eigenvalues. We have as a good approximation of 


= .01 (M) U (03.) U 


(6e) 


where 03 . is a diagonal matrix of the square roots of the eigenvalues of 
Eq. (6d). The matrix is given in Table 8. 
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The SOUSSA Program and the Fade Approximants 


The SOUSSA digital program (Refs. 8,9) computes generalized aero- 
dynamic forces for a wing of given planform executing sinusoidal oscillations 

in a fixed wing deflection mode shape. The generalized forces are computed 

cob 

at a given Mach number, m , for a given non-dimension frequency, k = , 

V 

a characteristic length, b , and a non-dimensional time variable T = — t . 

The generalized forces are given in terms of the force per unit dynamic pres- 
sure, q. Thus we have 



(m,k) + j 


Qj (m,k)^ 


W (j o>) 


( 7 ) 


where j = and W(j o>) is the non-dimensional deflection mode shape 

oscillating at the sinusoidal frequency, (U . 

hi the application we seek in this study, we desire to obtain the aero- 
dynamic forces in influence coefficient form. To obtain this, a pre-processor 
was developed by Prof. L. Morino to generate 17 unit impulse function modes for 
W(jw) corresponding to each of the 15 coordinates of our dynamical state w , 
plus the root section degrees-of-freedom, h(0) and a(0). Thus, for the ith unit 
impulse function mode we have imit deflection for the ith element and zeroes 
for the other fourteen elements. 


We obtain the ith column of the desired aerodynamic influence coefficient 
matrix from Eq. (7) 



(m,k) + j q^. (m,k) 


( 8 ) 


The SOUSSA coordinate system is positive deflection upward, and a 
positive angle of attack is trailing edge upward. Consequently, we have 




By choosing the characteristic length b to be unity, we have 



(9b) 


The SOUSSA generalized force is positive upward, and the generalized 
moment is positive for a positive Upward) force actii^ aft of the rotation 
axis. Thus, we have 





(9c) 


Finally, we have from the relation between time, t , and T 



The resulting unsteady lift, x (J oi)» for the non-dimensional sinusoidal 

P 

W (jw) vector is given by the Laplace transform of the last fifteen (15) rows 
of the matrix Eq. (6), 



- () f k i-F ) Vhj +j ikH^- (^k)%)4(Ja»} (11) 


and from Eq. 's (7) and (9c) we have 



(m,k) + j Qj 



( 12 ) 


It then follows that 

(j ^ k I-Fp)"^ (h^ + j Hg - (-^k)^Hg) = Qj^(m,k) + j Qj(m,k) (13) 
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To determine the desired constant matrices, independent of time and fre- 
quency, we have recourse to the Pade'^Approximants of Ref. 's 7,11,12,13. 

In what follows below, we lean heavily on the work of Edwards in Ref. 12. 

Since there are four (4) unknown matrices to be determined, we can, at most, 
satisfy only four (4) conditions. For the first condition, we choose to satisfy 
Eq. (13) at k = 0. We have 


Hi = -FpQu(m,0) 


(14) 


For the second condition, we choose to determine H from the real part of 

O 

Eq. (13) givii^ 

,2 N b 

H 3 = (-^) F (Qj^(m,k)-Qj^(m,0)J;+(— )Qj(m,k) (15) 


As k increases beyond bound, we have the piston theory limit (Ref. 14) 


= Q . X (in) (b = 1. ) 

3 V piston 


(16) 


The nonzero elements of Q . . , q.. , are given by (i = 1,7) 

^piston 13 




. 2 2' 
,4a / S - S+1 

m 


*^i, i+7 m 


, C.^ - a ^ - 2 X 

4a ( 1 1+1 o 



V7,i= \i+7 


q , q = j(l-Sx 

^1+7, ^1+7 m V 4 y ' o o' 


‘^,15 


3 3 

^ -S 

m 


(1 - X^) (1 + x^ - 2 x^) 


(16a) 
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^15,6 

^15,15 



1 - X. 




^ 6,15 


^13,15 

4 4 3 

. C_ - C„ 1 “ X. 

4g 6 7 1 

■ m 4 3 


(16a) 


where 


1.9431 


(X 

.8764016 -,3431794 


C. = wing chord in stream direction at the start of the ith section 

X = elastic axis in % of chord 
o 

(.4231) 

Xj^ = hinge line axis in % of chord 

(.80) 


The values of the are 


1 
1 

2 

3 

4 

5 

6 

7 

8 


C. (meters) 

.8353 

.7534 

.6784 

.6087 

.5390 

.4697 

.4205 

.3432 
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We choose as our third condition to match the imaginary part of Eq. (13) to 
the SOUSSA output for the flutter-reduced frequency, . Thus we have 



Fp Qj (m.k^) 


(17) 


Finally, we choose for a diagonal matrix to provide stable poles for the 

homogenous differential equation for the x variable. We choose for F 

P P 



when b, , b„ . . .b„ are the local semi-chords of the seven y. stations at 
12 7 1 

which the 15 coordinates are defined. The semi-chord used for the control 
surface (i = 15) is the semi-chord corresponding to the sixth (6) wing panel. 

In order to determine the best value of tr , a one dimensional search 
was undertaken to determine the open loop flutter analysis for the homogenous 
matrix differential equation (Eq. (6c) ) at a dynamic pressure of q = 5.36 kPa, 
a mach number of . 897, and a v of 136 meters/sec. The best value of 
proved to be a = +2. 249. 
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vn. Open Loop Flutter Analysis 


The results of the open loop flutter analysis of the DAST wii^ using 
the influence coefficient method was carried out by determining the eigen- 
values of the homogenous part of the differential matrix equation (6 c). 



r -N 

W ' 


I 

0 'X 


— ^ 

dt > 

w 

r 

= j -m"^k 

1 

1 

q m"^ \ { 

w 

1 


X 


h.-h„m"^d 

2 3 w 

F^ +qH„M"7 
p 3 / 

L vj 


(19) 


for different values of the dynamic pressure, q, for a fixed Mach number 
(m = . 897) and fixed airspeed ( v = 136 m/sec). 


The results of the study, shown in Figure 13, are to be compared to a 
similar plot taken from Ref. 1 shown in Figure 14. In order to illustrate a 
more detailed comparison of the pertinent flutter modes, we have Fig. 15 
which is a plot of the frequency and damping versus d 3 mamic pressure for the 
open loop system. 

The wind tunnel results as obtained from R. Doggett were 


m 

= 

.897 

q ^ 

= 

5. 36 kPa 

V 

= 

136 m/sec 

CC 

= 

8.0 Hz 


( 20 ) 


The comparison is shown to be good with both the wind tunnel data and 
the analytical prediction of Abel (Ref. 1). 
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Fig. 1 Wing Planform of DAST-ARWl 



TABLE 1 SPAR GEOMETRY DETAILS 


DISTANCE 
ALONG 
ELASTIC AXIS 
METERS 

A 

METERS 

B 

METERS 

C 

METERS 

t’ 

METERS 

0 




_ 

. 17145 

. 0331216 

. 0722122 

. 212344 

.00508 

.393192 

.021082 

. 059436 

.201422 

. 00508 

.78232 

.01905 

. 055372 

. 182372 

.00508 

1. 117346 

.017018 

. 050038 

.17526 

.00381 

1.452118 

. 015748 

. 044196 

. 150368 

.00381 

1.787398 

. 014097 

. 03937 

.128524 

.00381 

2. 122424 

. 013208 

. 03302 

.10287 

.00381 

2.426208 

.010922 

. 031242 

. 094488 

.00381 


B 



Fig. 2 Sketch of spar cross section along elastic axis 
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No. m. 
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4 1. 1049 

5 1. 35382 

6 1.60782 

7 1. 83642 

8 1.661414 
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Fig. 4 Coordinate System 



STRUCTURAL AXIS 




DEGREES OF FREEDOM AT SECTION Yj 


Fig. 5 Dynamical coordinate system 
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All grid points shown were restrained 
in DOF 3 and 5. Grid point 100 was 
also restrained in DOF 4. 


Figure 7. 


Point 111 is the location of the con- 
straint to the actuator torque. 

Points 97, 100 and 102 are constrained 
in DOF 3, 4 and 5. 


Grid Points Constrained in Iteduced SPAR Model 
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TABLE 3 

COMBINED BEAM AND PLATE ELEMENT MASS PROPERTIES 


Section 

No. 

Mass 

(kgm) 

Iy’ 

^ 2 
(kgm-m ) 

2 

(kgm^mi ) 

(kgm-in^) 

. % 2 
(kgm-m ) 

(kgm-m 

0 

7.1088 

.0262 

.0654 

.0913 

.0487 

.0429 

1 

2.5310 

. 0048 

. 0343 

.0382 

.0218 

.0173 

2 

1.8328 

.0026 

.0199 

.0220 

.0126 

,0100 

3 

1.3659 

.0017 

. 0128 

.0142 

'.0081 

.0064 

4 

1.1339 

.0012 

.0107 , 

.0116 

- .0066 

.0052 

5 

.9982 

.0010 

.0095 

.0103 

.0059 

. 0046 

6 

.7266 

. 0002 

.0058 

.0061 

.0034 

. 0026 

7 

.3341 

.0001 

.0023 

.0024 

.0014 

,0011 

* T" A 'T^TPTi's 


16.0313* 














* 

Compares with 15.8773 kgm calculated by SPAH 
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TABLE 4 


CONCENTRATED MASS PROPERTIES IN X-Y COORDINATE SYSTEM 


Section 

Mass 

Pitch 

Yaw 

Roll 

AX 

AY 

Pitch 

No. 

(kgm) 

Inertia 

2 

(kgm-m ) 

Inertia 

2 

(kgm-m ) 

Inertia 

2 

(kgm-m ) 

(m) 

(m) 

Inertia 

lY(grid point) 

, 2 
(kgm-m ) 

1 

1.6556 

.0734 

.0801 

.0134 

-. 0051 

+.0152 

.0734 

2 

1.0342 

. 0430 

.0443 

.0044 

+. 0147 

+. 0053 

.0432 

3 

1. 1975* 

. 0333 

.0344 

.0042 

+. 0239 

+.0036 

.0340 

4 

1.1612* 

.0267 

.0280 

.0034 

.0201 

.0028 

.0272 

5 

1.0161 

.0177 

.0192 

.0027 

.0218 

.0041 

.0182 

6 

1.4198 

.0200 

.0206 

.0047 

.0285 

.0000 

.0211 


.6350 

.0101 

.0127 

.0021 

-.0521 

-.0112 

.0118 


.5534 

.0056 

0 

.0056 

.2070 

-.0210 

.0293 


.4627 

0 

0 

0 

0 

0 

0 


9. 1355 

L , . ^ = I + m AX^ 

^(grid point) Y 

These values were altered to agree with SPAR computer model 

AX - Mass offset from grid point in X direction 
AY - Mass offset from grid point in Y direction 

These masses were added to cause flutter within the available wind tunnel dynamic pressure. 
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TABLE 5 MASS MATRIX 



42809?350E-03 





VIBRflTIONRL MODE. FREQ (HZ) .527873 XIO 


OJ 




o 

+ 




CJ 

LU^ 
Q_ ■ 
CO •' — ' 


- 31 - 


Fig. 9a Mode 1 (5 Hz) 
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SPEC WIND TUNNEL MODEL 



SPEC NIND TUNNEL MODEL 


VIBRATIONAL MODE, FREQ (HZl .601330 XIO 



SPEC WIND TUNNEL I10DEL 







TABLE 6 


COMPARISON OF VIBRATION MODES 



Vibration 

(Hz) 

Reduced 


Test 

SPAR 

Model 


Freq. 

666 DOF 

15 DOF 

1 

5.22 

5.279 

5.456 

2 

19.44 

18.90 

18.54 

3 

26.04 

26.01 

24.88 

4 

45.96 

44.34 

42.13 

5 

N/A 

60. 19 

61.67 

6 


73.87 

68.42 

7 


84. 96 

93.02 

8 


93.18 

105.45 

9 



126. 12 

10 



135. 01 

11 



146.70 

12 



204.40 

13 



297.65 

14 



420. 15 

15 



894. 27 
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Fig. 13. Open Loop Flutter Analysis of the DAST wing for fixed 
Mach Number = . 897 and v = 136. m/sec 
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Figure 14. Dynamic-pressure root locus at M = 0.90 (system off) . 
Arrows Indicate Increasing dynamic pressure. 
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